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Abstract 

Strongly correlated quantum impurity problems appear in a wide variety of contexts rang- 
ing from nanoscience and surface physics to material science and the theory of strongly 
correlated lattice models, where they appear as auxiliary systems within dynamical mean- 
field theory. Accurate and unbiased solutions must usually be obtained numerically, and 
continuous-time quantum Monte Carlo algorithms, a family of algorithms based on the 
stochastic sampling of partition function expansions, perform well for such systems. With 
the present paper we provide an efficient and generic implementation of the hybridization 
expansion quantum impurity solver, based on the segment representation. We provide 
a complete implementation featuring most of the recently developed extensions and op- 
timizations. Our implementation allows one to treat retarded interactions and provides 
generalized measurement routines based on improved estimators for the self-energy and 
for vertex functions. The solver is embedded in the ALPS-DMFT application package. 
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Classification: 7.3 
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External routines/libraries: ALPS [UEIE], BLAS [HE], LAPACK [6], HDF5 [7J . 
Nature of problem: Quantum impurity models were originally introduced to describe a magnetic 
transition metal ion in a non-magnetic host metal. They are widely used today. In nanoscience 
they serve as representations of quantum dots and molecular conductors. In condensed matter 
physics, they are playing an increasingly important role in the description of strongly correlated 
electron materials, where the complicated many-body problem is mapped onto an auxiliary 
quantum impurity model in the context of dynamical mean-field theory, and its cluster and 
diagrammatic extensions. They still constitutes a non-trivial many-body problem, which takes 
into account the (possibly retarded) interaction between electrons occupying the impurity site. 
Electrons are allowed to dynamically hop on and off the impurity site, which is described by a 
time-dependent hybridization function. 

Solution method: The quantum impurity model is solved using a continuous-time quantum 
Monte Carlo algorithm which is based on a perturbation expansion of the partition function in 
the impurity-bath hybridization. Monte Carlo configurations are represented as segments on 
the imaginary time interval and individual terms correspond to Feynman diagrams which are 
stochastically sampled to all orders using a Metropolis algorithm. For a detailed review on the 
method, we refer the reader to [8]. 

1. Introduction 

Quantum impurity models describe a set of correlated sites or orbitals embedded 
in a bath of non-interacting states. Quantum impurity models appear in a range of 
contexts, including magnetic impurities embedded in a non-magnetic host material [5], 
nanoscience, where they are used to describe quantum dots and molecular conductors 
0, and surface science [TU], for the description of molecules adsorbed on a substrate. 
Quantum impurity solvers are also an essential ingredient of the dynamical mean -field 
(DMFT) [n] [12l US] HI] approximation to correlated lattice systems, which has had 
enormous success in recent years in the simulation of correlated material systems |15[ I16[ 
[17] and lattice models [TB] . 

With this article we provide a description and a state-of-the-art implementation of 
the continuous-time [1911201, [2"T] 'hybridization expansion' quantum Monte Carlo impurity 
solver for density-density interactions |22j . Our implementation includes in particular 
the important numerical and conceptual advances developed over the last few years: 
improved estimators |23j . frequency and Legendre measurements |24j . measurement of 
vertex functions [3T] , treatment of retarded interactions [25] [2B] , and parallelization to a 
large number of cores [T] [2] [3] [27J . 

General fermionic impurity models have the form -ff; mp = H\ oc + -ffbath + ^hyb > where 
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The term ffi oc corresponds to the impurity with level energies and intra-orbital hoppings 
described by E and the interaction terms parametrized by U (roman indices label the 
different interacting orbitals including spin), -ffbath describes the non-interacting bath 
with quantum numbers k and spin/orbital index a. The hybridization term -ffhyb rep- 
resents the exchange of electrons between the impurity and the bath, parametrized by 
the hybridization matrix V£* b . All the relevant properties of the bath are encoded in the 
hybridization functions 

A ab (tw„) = k _ ■ (4) 

In a DMFT calculation, these parameters are determined self-consistently. 

The observable of primary interest is the impurity Green's function G, defined a^] 

G ab (r-r'):=-(T T c a (r)cl(r')). (5) 

The Green's function contains the information about all single-particle properties of the 
model, including the spectral function and the particle density. Higher-order correlation 
functions are needed for the measurement of susceptibilities and the calculation of the 
(reducible or irreducible) vertex. 



2. Hybridization Expansion 

In recent years, a new class of efficient Monte Carlo techniques has been developed 
for solving quantum impurity models: the continuous-time impurity solvers [2TJ [L9\ [2Ql 
[221 HS1 [2H1 [50]. The hybridization expansion [221 HE] approach that we describe here 
is particularly well suited for the single- and multi-orbital impurity problems that typi- 
cally appear in single-site DMFT calculations. In this approach, the partition function 
Z = Tr £ ;, c [e~^- HhTlp ] is expanded in powers of the hybridization term -ffhyb- Monte Carlo 
configurations consist of a sequence of creation and annihilation operators, and thus 
represent a sequence of hybridization events (electrons hopping from the bath to the im- 
purity or back into the bath). The Monte Carlo weight w({Ti}) of such a configuration 
is the product of two factors, 

w ({ T i}) = moc({Ti})w hyh {{Ti}), (6) 

where w\ oc ({Ti}) is the trace over the impurity states for the specific sequence of impurity 
creation and annihilation operators, and ui n yb ({"?!}) is the trace over the bath states. 
Because the bath is non-interacting, the latter can be expressed as the determinant of a 
matrix M , whose elements are hybridization functions connecting the different pairs 
of creation and annihilation operators. 

The local contribution is evaluated explicitly. For a generic model, i. e. if the Hamilto- 
nian is not diagonal in the occupation number basis, this is a computationally expensive 
procedure which scales exponentially with system size [28, 29, 31j, even though substan- 
tial progress has been made to make these simulations affordable [3H [5TJ [33J . An 



iWe adopt here the 'many-body' definition of G with G(0-|_) < as opposed to the 'Monte Carlo' 
definition where G(0+) > 0. 
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important simplification occurs if the interaction term is restricted to density-density 
interactions. In this case, the atomic part of the Hamiltonian takes the form 

Hioc = Ean - + J2 uabn * n b (7) 

a ab 

and the local trace can be evaluated in polynomial time [22j . 

3. Implementation 

We aim to provide a fast but general implementation of this algorithm. We include 
a formalism for retarded interactions, efficient measurements of single-particle Green's 
functions in imaginary time in both the Matsubara and Legendre basis, measurements 
of two-particle Green's functions, improved estimators for the self-energy and vertex 
functions, and energy and sector statistics measurements, as well as density-density cor- 
relation functions. We briefly describe each of these features in the following. More 
detailed explanations can be found in a recent review |21j and the original papers. 

3.1. Retarded Interactions 

Retarded interactions appear in models with electron-phonon coupling (e.g. the 
Hubbard-Holstein model |34| ) or in realistic material simulations with dynamically screened 
Coulomb interaction. They also arise in the context of extended dynamical mean-field 
theory [351 ES] an d its diagrammatic extension, the dual boson approach [37]. An efficient 
technique to treat retardation effects in the hybridization expansion has been described in 
Refs. [25] and [26] . In order to treat a model with a frequency-dependent local Coulomb 
interaction U(u>) one has to (i) set the instantaneous interactions to the screened value 
U(oj = 0) and (ii) introduce an effective retarded interaction K(t) between all pairs of 
creation and annihilation operators (irrespective of flavor). The explicit expression for 
the function K(t) in the interval < r < (3 is [26] 

K(r)= r*±™M [WuiT )-W um , (8) 
Jo * w 

with W u (t) = cosh [(t - § )u] / sinh ["f] . The sign of the retarded interaction depends 
on the type of operators (creators/annihilators) which are connected, so the additional 
weight factor of the Monte Carlo configuration becomes 

Screen (M) = >>, (9) 

where the sum is over all the operator pairs i, j and Si = +1 if operator i is a creation 
operator and —1 if it is an annihilation operator. 

Frequency-dependent (retarded) interactions can thus be incorporated very easily into 
the hybridization expansion method. If a new pair of creation and annihilation operators 
(segment) is introduced into a Monte Carlo configuration, one computes the change in 
wioc (using U(lu = 0)) and the change in w scrccn . For the latter, one has to evaluate 
the retarded interactions between all the operator pairs which involve at least one of the 
newly inserted creation and annihilation operators. The computational effort for a local 
update is O(n) and therefore negligible compared to the evaluation of the determinant 
ratio, which is 0(n 2 ). 
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3. 2. Green's function measurements 

3.2.1. Imaginary time measurement 

For a given inverse hybridization matrix M, the Green's function, Eq. |5]), is measured 
on the interval [0, 0] as 

G ab (T) =(^j Y.M lja 5-{T,T e a -Tl)8 aa 8 b }j (10) 

where (. . .) denotes the Monte Carlo average and S~(t, t') := sgn(r')5(r — r' — 6(—t')/3) 
accounts for the antiperiodicity of G(r). The r e (t s ) denote the times of annihilation 
(creation) operators. The Green's function is measured on a discrete grid representing 
the continuous variable r. Since the complexity of the algorithm does not depend on the 
grid resolution, this grid can be chosen arbitrarily fine. 

3. 2. 2. Matsubara frequency measurement 

It is often convenient to use the Matsubara frequency representation of the Green's 



function. Carrying out the Fourier transformation of Eq. ( 10 ) analytically, the Mat- 
subara coefficients G(iu n ) of G(t) are measured directly, without the need of a time- 
discretization grid: 



G ab {iu n ) =( — VM^c»'K- T ?)uJ ■ (11) 




3.2.3. Orthogonal polynomial representation 

A measurement procedure based on an expansion of the Green's function on the in- 
terval [0,(3] in terms of orthogonal polynomials was introduced in Ref. [23]. We focus here 



on Legendre polynomials. Again, the transformation of the measurement rule Eq. (10) 
is carried out analytically. The Legendre transformation is unitary and can be written 
in terms of a matrix multiplication G(iu n ) — ^2i>oT n iGi, with expansion coefficients 




Gata = { ^M^P^tI - T$)5 aa 5 bf) ) (12) 

where P/(r) := Pi[x(t)] for r > and — Pi [x(t + 0)\ for r < and x(t) = 2t//3 — 1 maps 
the interval [0,/3] to [—1, 1]. Observable estimates exhibit a clear plateau as a function 
of I, making a basis truncation well controlled. The coefficients G/ generally decay faster 
than any power of l/l, leading to a highly compact representation of observables. As the 
noise is mostly carried by high-order coefficients, this representation acts as an efficient 
noise filter. 



3.3. Improved estimators 

An important physical quantity is the impurity model self-energy. It is usually com- 
puted by inverting the interacting and non-interacting Green's functions of the impurity 
according to the Dyson equation 

S(icj n ) = Gq 1 {iu n ) - G _1 (zw„). (13) 
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This inversion amplifies statistical noise, in particular at high frequencies where the 
difference between the Green's functions G and Go (which both decay as l/u>„ with 
constant errors) is small. 

This problem has been resolved recently |23) . Using the equation of motion, the 
product of the self-energy and Green's function can be expressed in terms of an additional 
correlator F^t-t 1 ) := -(T T c a (T)cl(r , )n j (T / )) (so-called "Bulla's trick" [38]): 

(GS) ni K) = \ YJPih + U bj )Fl b {iuj n ). (14) 

3 

The additional quantity can be computed from G at minimal extra cost, both in the 



Matsubara and Legendre basis (cf. Eqs. (Ill, (12)). Combining the Legendre filter 
with the improved estimator yields self-energies with unprecedented accuracy (see Sec. 
[4] for an illustration). The current implementation has been generalized to evaluate the 
improved estimators in the presence of retarded interactions. Details will be published 
elsewhere. 

3-4- Two-particle Green's functions 

Two-particle Green's functions and the related vertex function have gained impor- 
tance in recent years, partly due to an increase in computational resources but also 
through the development of new methods. They serve as input for the calculation 
of susceptibilities within DMFT and provide the basis for novel extensions of DMFT 
[33 HOI EH- To meet these requirements, we provide measurements for the two-particle 
Green's function 

G {2) (T a ,T b ,T c ,T d ) := (T T c a (T a )cl(n)c c (T c )c d {T d )) . 

These functions are measured in frequency space, as a function of two fcrmionic 
frequencies and one bosonic frequency. We further provide the measurement of the 
correlation function 

H 3 {T a ,T bl T c ,T d ) := 

(r r n J (r a )c a (T Q )c] ) (T f) )c c (T c )c^(r (i )). (15) 
which allows one to accurately compute the vertex function from an improved estimator 



expression in analogy to Eq. (14). For details, see Ref. 



3.5. LDA+ DMFT interface 

The LDA+DMFT method [17l[T6] provides an interface between band structure meth- 
ods and many-body theory: A local density approximation (LDA) band structure is ob- 
tained from a band theory calculation, to which two-body correlations (typically static 
Coulomb repulsion and Hunds coupling terms) are added. Our code has been designed 
specifically with these applications in mind, and it provides an interface to specify in- 
teraction matrices and double-counting corrections. In the context of LDA+DMFT, 
off-diagonal hybridizations and non-density-density interactions may become important. 
While related algorithms exist to treat them [2H 131] > and substantial progress has been 
made to make these simulations affordable [25] |32J [5TJ [33] , the present implementation 
is restricted to diagonal hybridization functions and density-density interactions. 
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3. 6. Sector statistics 

A histogram of the atomic states occupied in the course of the simulation may provide 
important physical insights |29j . An impurity with a single orbital can assume one of 
four different states: empty (|0)), singly occupied (|f) or \\)) or completely filled (|t4-))- 
The algorithm measures the fraction of time the impurity spends in any given atomic 
state and collects a histogram. A typical example is shown in Fig. [6] 

3. 7. Parallelization 

Continuous-time Monte Carlo codes are straightforwardly parallelizable by running 
a different random walk on each core. The cost of thermalization and the non-parallel 
sections of the code are negligible, therefore efficient implementations scale up to a com- 
paratively large number of CPUs [27 . The present implementation of the hybridization 
expansion code supports MPI parallelization using the next-generation ALPS scheduler 
PP. Efficient use of collective communications and state-of-the-art binary data storage 
[7] mean that for typical applications the code scales to more than 1000 cores without 
noticeable overhead. 

3.8. Python interface 

The solver can be built either as a standalone executable or as a Python module. 
Using the solver as a Python module allows one to perform all computationally expensive 
operations in C++, while other aspects of the calculation (e.g. a DMFT self-consistency 
cycle, data evaluation, or application specific computations) can be written in Python. 
All data storage is implemented using the popular hdf 5 library, with an interface both 
to C++ and Python. 

4. Examples 

Single- orbital Anderson impurity model 

As a first example, we consider the single-orbital Anderson impurity model. In the 
following, we always use a semielliptical density of states with bandwidth equal to At, 
and set t = 1. Fig. [T] shows the perturbation-order histograms for U — 2 and different 
temperatures. The mean perturbation order shifts to larger values with decreasing tem- 
perature. Note that the average perturbation order yields the kinetic energy divided by 
temperature [29j . 

As an example for a non-trivial observable, we measure the local spin-spin suscepti- 
bility 

rf> 

Xdd(iu m ) ■= / dT{S z (T)S z (0))e WmT (16) 
Jo 

directly in frequency. Fig [2] shows its static component Xdd{w — 0) times temperature 
(i.e., the effective local moment) versus temperature. As the temperature is lowered a 
local moment is formed, the magnitude of which increases with U. As the temperature is 
lowered further, the moment decreases as a result of screening due to the Kondo effect. 
The screening occurs on a scale set by the Kondo temperature Tr-. With the present 
solver, it is possible to reach very low temperatures. 
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Figure 1: (Color online) Perturbation-order histograms for the single-orbital Anderson model with 
U = 2 for different temperatures. 

4-2. Retarded interactions 

As a second application, we consider a single-orbital model with a retarded inter- 
action, corresponding to a Holstein-Hubbard model ("plasmon" frequency electron- 
boson coupling strength A), which we solve within DMFT on the Bethe lattice. The 
retarded interaction function K(t) and its derivative are given by 

K(t) = -4{ cosh M/V2 - r)]/sinh(w /3/2)} + c, 

\2 

K'(t) = +— sinhM^/2 - T)]/sinh(w /3/2), 

LOQ 

where the constant c is chosen such that K(0 + ) = 0. Fig. [3] shows the electronic 
self-energy for fixed bare interaction U = 8i and fixed screened interaction U scr — 
U — 2A 2 /ujq = 'it for different screening frequencies u>o at temperature T = 0.02i. For 
large screening frequencies the system is metallic, while it undergoes a metal-insulator 
transition as ujq is lowered. The corresponding results for the spectral function of this 
model are shown in Ref. |26| . 
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Figure 2: (Color online) Effective local moment of the single-orbital Anderson model as a function 
of temperature. As the temperature is lowered the local moment is formed and then decreases due to 
Kondo screening. 



4.3. DMFT for Multi-Orbital Models 

The code is well suited for studying multi-orbital problems in the density-density 
approximation. As an example we consider a two-orbital model with interaction 

22 U ab n a n b = U ^ n at n al + U '*^2 n h<r n 2,-T 
ab a— 1.2 a 

+ (£/'- J)Y,m,*ri2,„- (17) 

Here a and a denote spin and orbital indices, U and U' = U — 2 J are the intra- and inter- 
orbital Coulomb interaction parameters and J is the Hund's rule coupling coefficient. 
We solve the model on the Bethe lattice with bandwidth At at temperature T/t = 1/50, 
U/t = 8 and relatively large Hund's coupling J/U = 1/6. We consider only paramagnetic 
and paraorbital solutions. This model has been shown to exhibit a spin-freezing transition 
as a function of filling [23] . This manifests itself for example in the spin-spin correlation 
function (see Fig. |4|, which is small for large times (i.e. near j3/2) in the Fermi liquid 
phase, but approaches a large non-zero value in the frozen moment phase. 

As a second example we show the imaginary part of the self-energy for different fillings 
across the transition (Fig. [5]), determined using the improved estimator measured in two 
different bases: the Matsubara basis and the Legendre basis (Ni = 80 coefficients). 

9 
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Figure 3: (Color online) Self-energy of the Hubbard Holstcin model with (7 = 8, fixed screened 
interaction U BC r = 3, temperature T = 0.02 and different phonon frequencies uiq. The self-energy has 
been measured using improved estimators in the Legendre basis. 

The imaginary part of the high-frequency tails Stall (iw n ) = + S^/(«w„) is shown 
for comparison. It has been calculated from the orbital densities (n*) and equal-time 
density-density correlation functions (njjij) using the expression 

=^U ik U il {{n k n l ) - (n k )(m)). (18) 

kl 

One can see that the Legendre representation efficiently filters the Monte Carlo noise 
and correctly reproduces the high-frequency tail. 

For the same parameters, we show in Fig. [6] the sector statistics, which measures 
the fraction of time the impurity spends in a given atomic state. One can see that 
for the lowest filling, the impurity essentially is either completely empty ( |0000) ) or 
singly occupied (|1000), |0100), |0010) or |0001)). Doubly or higher occupied states 
are unlikely, so correlation effects arc weak. Indeed, at these parameters the system 
has previously been identified as a weakly correlated Fermi liquid |23| . For a higher 
filling, (n) = 1.29, the system is on the boundary between the Fermi liquid and the 
frozen moment phase. While the probability of the singly occupied states remains almost 
unchanged, the contribution of the empty state becomes small, while the doubly occupied 
high-spin states (| 1010) and 1 0101) ) become significantly populated. This indicates the 
formation of the spin 5 = 1 moment. The low-spin doubly occupied states (| 1100) and 
|0011)) are suppressed due to the large intra-orbital Coulomb repulsion U, while |0110) 
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Figure 4: (Color online) Spin-spin correlation function as a function of imaginary time for various values 
of filling. Two-orbital model, U/t = 8, J/U = 1/6, fit = 50. 

and 1 1001) are suppressed due to U' . The latter have higher weight since U' < U . For 
even higher filling ((n)=1.53), deep in the frozen moment phase, the high-spin states 
gain further importance at the expense of the singly occupied states. Additional results 
for this model can be found in Ref. [23"] . 

5. Installation and Usage 

5.1. Installation 

The program is provided as part of the ALPS [lj package and can be downloaded from 
|http://alps .comp-phys.org. Release milestones and binary nightly builds as well as the 
source code are available as a package and via anonymous svn access. The code binary 
hybridization can be found in the ALPS binary directory alps_root/ bin. Detailed 
documentation with an overview of all parameters and example scripts is located in 
alps_root/doc as hybdoc.pdf . 

5.2. Usage 

The code is delivered together with tutorial examples which illustrate the usage of 
the program and allow one to reproduce the data presented here. For a quick inspection 
and test of the program, do the following. After installation of the library 
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Figure 5: (Color online) Self-energy of a two-orbital model (U/t = 8, J/U = 1/6, T/t = 0.02) for 
different fillings, calculated using the improved estimators. Dots represent data from the Matsubara 
measurement, whereas the result from the Legendre measurement is shown by thick solid lines. Thin 
solid lines correspond to the high-frequency tails. 



• From the directory alps_root/tutorials enter directory hybridization-01-pyth.on. 

• Inspect the file tutoriall .py, which describes the parameters. 

• Run alps_root/bin/alpspython tutoriall .py. This run will take about 30 sec- 
onds. 

• Plot Gt.dat and compare it to Gt.dat in the directory exact_diagonalization, 
which contains exact results. 

• Plot the file orders.dat, which contains the expansion orders. 
For a more extensive look at the program, do the following. 

• Enter directory hybridization-02-kondo. 

• Inspect the file tutorial2.py, which describes the parameters. 

• Run alps_root/bin/alpspython tutorial2.py. The script will start several con- 
secutive runs which each take a few seconds. After completion, a plot appears. 
Compare this plot to Fig. [2j 

Two more tutorials generate the data in Figs. [3][6j The extensive documentation in the 
subdirectory doc provides more information and additional examples. 
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Figure 6: (Color online) Sector statistics for the two-orbital model (U/t = 8, J/U = 1/6, T/t = 0.02) 
showing the probability of finding the impurity in a given atomic state. The states are labeled by the 
occupation numbers of the four distinct spin-orbital states. 



6. Conclusions 

In conclusion, wc provide here a state-of-the-art implementation of the hybridization 
expansion CT-QMC method which takes into account recent advances in methodology. 
Our code package also includes a set of examples that provide a pedagogical introduction, 
as well as Python scripts to reproduce the examples in this paper, and allow the user 
to test the implementation. The code has been written specifically with users from the 
LDA+DMFT and nano-science community in mind, with the aim of providing both a 
starting point for students and a high-performance implementation that can be used in 
a production environment. 
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